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Abstract 



Motivated by recent experimental studies of Bodenschatz et al. [E. Bodenschatz, J.R. 
de Bruyn, G. Ahlers and D.S. Cannell, Phys. Rev. Lett. 67, 3078 (1991) ], we present 
a numerical study of a generalized two dimensional Swift-Hohenberg equation to 
model pattern formation in Rayleigh-Benard convection in a non-Boussinesq fluid. It 
is shown that many of the features observed in these experiments can be reproduced by 
this generalized model that explicitly includes non-Boussinesq and mean flow effects. 
The spontaneous formation of hexagons, rolls, and a rotating spiral pattern is studied, 
as well as the transitions and competition among them. Mean flow, non-Boussinesq 
effects, the geometric shape of the lateral wall, and sidewall forcing are all shown to be 
crucial in the formation of the rotating spirals. We also study nucleation and growth 
of hexagonal patterns and flnd that the front velocity in this two dimensional model is 
consistent with the prediction of marginal stability theory for one dimensional fronts. 



2 



1 Introduction 



One of the most striking examples of spatio-temporal self-organized phenomena in 
nonequilibrium systems is the rotating spiral states observed in chemical and bio- 
logical systems [Q]. It is remarkable that such time dependent but macroscopically 
coherent states can be sustained in systems that are not in equilibrium. The Belousov- 
Zhabotinsky reaction 0, for example, has received considerable attention as an ex- 
ample of a chemical wave propagation. Spiral patterns in this system result from the 
couphng of reaction and transport processes. 

Recently, similar rotating spiral patterns have been observed in Rayleigh-Benard 
convection in non-Boussinesq fluids in large aspect ratio systems 0]. According to 
the classical work of Busse [Q, the first bifurcation from the conducting state is 
to a convective state of hexagonal symmetry. Convective cells form a stationary 
honeycomb structure. Further away from threshold, the system undergoes a new 
bifurcation to a state comprising parallel convective rolls (roll patterns are the only 
patterns predicted and observed within the Boussinesq approximation. The existence 
of of a stationary pattern of hexagonal symmetry is a direct consequence of deviations 
from the Boussinesq approximation). The predicted bifurcation from hexagons to 
rolls is direct, so that the fluid is expected to evolve to a stationary patterns of rolls. 
Recently, however, experiments on convection in CO2 gas by Bodenschatz et al. |0 
show that the system has a tendency to spontaneously form rotating spirals. These 
rotating states are long lived, and do not decay to the expected pattern of concentric 
rings (in a circular geometry). Furthermore, depending on the value of the Rayleigh 
number, spirals with a different number of arms have been observed. 

Rayleigh-Benard convection in monocomponent fluids is governed by the full three 
dimensional fluid equations. Because of the difficulty in solving the three dimensional 
initial value problem posed by the fluid equations, we and others have focused on the 
study of simpler two dimensional model equations. An example of such models is 
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the so-called "Swift-Holienberg" (SH) equation, which is asymptotically equivalent 
to the long distance and long time behavior of the fluid equations near onset of 
convection and in the Boussinesq approximation 1^, A great deal of theoretical 

and numerical work on this latter model has been done by Cross , and by Greenside 
et al. P], but much remains to be done in non-Boussinesq systems, even within the 
framework of a SH-type equation. 

We study in this paper a generalized SH equation that includes both a quadratic 
nonlinearity and coupling to mean flow effects. The values of the parameters that en- 
ter the equation have been chosen to be in the range appropriate for the experiments 
of Bodenschatz et al. We find that stable rotating spirals are spontaneously formed 
during the hexagon to roll transition, in agreement with the experimental observa- 
tions. The quadratic nonlinearity in the equation is responsible for the rotational 
symmetry breaking and leads, by itself, to stationary spiral patterns. When coupling 
to mean flow in included, and therefore when the model equation is not potential, the 
same transition leads to rotating spirals instead. However, spirals are not obtained 
if one sets the quadratic term equal to zero but keeps the mean flow term. Sidewall 
forcing also seems to be essential in obtaining this pattern. Otherwise, rolls that 
are locally perpendicular to the sidewall appear and no uniformly rotating state is 
observed. Finally, once the spiral state is formed, it is unstable to the removal of the 
quadratic nonlinearity in the equation. The spiral state quickly decays to a set of 
concentric rings. 



We note that Bestehorn et al. |10| have reported a numerical study of this very 
same model to study the rotation of a spiral pattern, but their study is limited to the 
special case in which the initial configuration is a spiral. The fundamental question 
addressed in our work is, in contrast, how the spiral pattern is spontaneously formed 
during the hexagon to roll transition. 

In Section II, we give a brief description of the theoretical model used to describe 
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pattern formation in non-Boussinesq systems, and in Section III we provide a detailed 
description of the numerical method for solving the model equation. In Section IV, 
we present various numerical results, show in detail the hexagonal and spiral patterns 
and compare our results with the experiments of Bodenschatz et al. In Section V, we 
present a brief summary. 



2 A two dimensional model equation for convec- 
tion in non-Boussinesq fluids 

A great deal of our current understanding of the spatial and temporal properties 
of convective patterns near onset has been obtained by using model equations such 
as the Swift-Hohenberg (SH) equation, which in dimensionless variables takes the 
following simple form |lT], 

dijj{f, t) 



(1) 



dt 

where r is a two dimensional vector lying in the (x, y) plane (the convective cell is 
parallel to this plane), and e is the control parameter. This equation describes the 
evolution of a single scalar field ip which is commensurate with the convective rolls. 
When e > 0, this equation has roll-like solutions with a wavenumber q = 1. This 
equation can also be written in potential form, 

dt~ 5r ^ ' 

where £ is a Lyapunov functional. The dynamical evolution of is naturally inter- 
preted in terms of the minimization of this functional. Therefore a potential system 
cannot display either the oscillatory or aperiodic time dependence observed in the 
experiments. Although potential systems are quite important for the understanding 
of the features of the emerging patterns, an interesting question one may ask is which 
flow component is necessary to describe the observed non-potential behavior. One 



such component is the so-called large scale mean flow, which was first proposed by 
Siggia and Zippelius [|T2 |. 



Large scale mean flows are composed of one particular spatial harmonic of the basic 
roll pattern. In the case of free-free boundary conditions, this harmonic corresponds 
to a flow that is independent of 2;, the coordinate normal to the convective cell plates. 
This flow is not damped in the (unrealistic) case of free-free boundary conditions, and 
is only slightly damped in the more realistic case of rigid-rigid boundary conditions. 
Since the nonlinear term in the fluid equations that gives rise to large scale mean 
flow is not V ■ V6, but v ■ Vw, the magnitude of the large scale flow is inversely 
proportional to the Prandtl number. Moreover, for a perfect straight roll pattern, no 
large scale mean flow is generated; it appears only when rolls are bent or modulated. 
In that case, the characteristic length scale of these flows is large compared with 
the roll wavelength. For example, the large scale mean flow contribution to the SH 
equation jl^ , |l3|| has been shown to play a key role in the onset of weak turbulence 



in Boussinesq systems [14 



We next describe a two dimensional model of convection in a non-Boussinesq fluid. 
We use the two dimensional generalized Swift-Hohenberg (GSH) equation ^, |T3 



given by Eqs. (|^) and (H) below, which we solve by numerical integration. Our model 
in dimensionless units is deflned by, 

d'ilj{f, t) 



dt 



e- V^ + 1 



d_ 

di 



where, 



with boundary conditions. 



(3) 
(4) 

(5) 



i\B = n- VeiB = V^ls = n ■ VV^Ib = 0, 



(6) 
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where h is the unit normal to the boundary of the domain of integration, B. If the 
couphng to large scale flow is dropped {gm = 0), Eq. is potential. Also this equa- 
tion reduces to the dimensionless SH equation (Eq. ([l|)) when the coupling coefficients 
92 = Qm = 0. As is the case in Eq. (|^), the scalar field ip{f,t) is proportional to a 
linear combination of the fiuid temperature modulation and vertical velocity fields at 
a point r in the midplane of the fiuid layer, parallel to the upper and lower walls of 
the convective cell, ^(r, t) is the vertical vorticity potential |jl2|, The quantity 

e is the control parameter. Pr is the Prandtl number of the fiuid, and iy are two 
unit vectors parallel to the x and y direction respectively, and is an unknown con- 
stant. A phenomenological forcing field / has been included in Eq. to simulate 
lateral sidewall forcing produced by horizontal temperature gradients present in the 
experiments. As in earlier studies |16|, we have varied the strength and spatial 
extent of / in order to match the experimental observations. In order to estimate the 
values of the various dimensionless parameters that enter the GSH equation in terms 
of experimentally measurable quantities, we have derived a three mode amplitude 
equation (see the Appendix for details). From the experiments described in reference 
IQI, we estimate that g2 = 0.35, gm = 50, = 10 and Pr = 1. The value of e used in 
the numerical calculation is related to the experimental value eexpt by eexpt = 0.3594e. 



3 Numerical method 



The numerical method for solving the GSH equation (Eqs. (^ - (|§)) is based on the 
elegant work by Greenside et al. 0, and Bj0rstad et al. [T^. We sketch below their 
numerical scheme. The key step is to recognize that the GSH equation can be solved 
by the repeated solution of the linear constant coefficient biharmonic equation for a 
function u, 

(V'V' + aV' + b)u = /i, (7) 



7 



with boundary conditions, 

m|r = /2, n-Vu\R = f3, (8) 

on a circular domain R. The constants a and b are real numbers, n ■ V denotes the 
normal derivative taken at the boundary of the domain, and /i , /2 and /s are given 
functions. For the case of rigid boundaries, we have /2 = /a = 0. 



Consider an implicit backward Euler discretization scheme in time, to yield the 
following finite difference set of equations for Eqs. (0)-®, 

+ + + Ar) . V.0(r + Ar) = 

L^{r + Ar) - g^ij^r + Ar) - ij^{r + Ar), (9) 



+ Ar) - e(r)] 



C2 



e(r + Ar) 



Ar 

- Wipir + At) X V(V2^(r + Ar)] ■ e^, (10) 

where iP{t) and ^(r) are the known solutions of Eqs. ®-(§) at time r, Ar is the 
time step, L = e — (1 + V^)^ is a linear biharmonic operator, and iIj{t + Ar) and 
^(r + Ar) are the unknown implicit solutions at the next time step (r + Ar) (we have 
temporarily suppressed the spatial arguments). We solve Eqs. (^) and (|10D by using 
a multi-iteration Gauss-Seidel scheme. We first assume that iIj{t + At) and ^(r + Ar) 
are obtained by successive approximations of the form, 

V^(r + Ar) ~ Vfc+i = ^fc + 4, withV'o = ^/'(r), (11) 

e(r + Ar) ~ ^fc+i = a + ^^fe, witheo = e(r), (12) 
where ipk and C,k are the approximations at the k — th iteration. The equation satisfied 
by the so-called outer correction fields, 6k and 9k, can be obtained (assuming || 6k \\ 



<^ II ipk II and II 6k \\ <^ \\ ^k \\ in the maximum norm). By substituting Eqs. ([TT]) 
and (0) into Eqs. (^ and ([ToD , and linearizing them with respect to 6k and 6k, we 
then obtain the standard Gauss-Seidel iteration scheme for the unknown corrections 
6k and 6k, 



6k = ^f^^ + QmUk ■ V^fc - (L^fc -i^l- g^i^l), (13) 



PrAr 



(14) 



with h"^ 



PrAr' 



The right hand sides of Eqs. (|T^) and (0) are the k — th outer 
residuals, r'^ter{k) and rl^terik) which measure the extent to which Eqs. (Q) and (|T0|) 
are satisfied by the k — th order approximation. Given the residuals, we solve for the 
outer corrections, 6k and 9k, and then obtain a better approximation to ?/'(r + At) 
and ,^(r + Ar). Iteration continues over the index k, until both the outer residuals 
and the outer corrections are small compared to ipk and ^k, that is. 



max(|| 4 II, II rtuter{k) \\) < e^e/dl i'ik) \\) + Cafes, 



(15) 



and, 



maxfll 9k 



rluterik) \\)<eU\\ m ID + eafes, 



(16) 



where e^ei and eahs are the relative and absolute error tolerances, chosen to be 0.1 and 
10"'^ respectively in our calculations. When the convergence criteria are satisfied, ipk+i 
and ^k+i are set to be ip{T + Ar) and ^(t + Ar). In the numerical solution, we have 
found that Ar > Ar^ax — 1-8 will cause a numerical instability, so we have chosen 
Ar < ATmax- (Ar ^1.0 during the initial transient at onset). For a given ipk+ii 



Eq. (|14|) has the exact form of Eq. (|^, with 



0, /i 



\k), 



p.f./\T^ ^ "-'^ Ji 'outer 

/2 = /s = 0, and can be solved rapidly and accurately. Eq. (14) is almost of the 
form of Eq. (|^ except for the non-constant term —2g2ipk — the left hand side 

operator. We can reduce this equation to the desired constant coefficient biharmonic 
form by assuming a successive approximation of the form: 



6k — 6k,m+l — 6k,m + Vn 



(17) 



where the inner correction field, rim{x,y), is assumed to be small compared to 6k,m, 
the m — th approximation to 6k- By substituting Eq. ([T7|) into Eq. (^) and solving 
for rjm by approximating the non-constant term acting on 1]^ as a constant C, we 
obtain 



L- — + C 
At 



■4) 

^ outer 



{k)- 



lAT 



The right hand side is the m — th inner residual, rf^j^^^{k,m) of Eq. (|T8|). This 
measures the extent to which Eq. (|13]) is satisfied after m iterations. The constant 
C = —2(72 < ■i/'fc > — 3 < ^/^^ >, where <> denotes a spatial average over the entire 
system. Now Eq. (p!8|) has the exact form of the constant coefficient biharmonic 
equation given in Eq.(0), with a = 2, b=l + -^ — e — C, fi = — r^„gr(A;, m), and 
/2 = /s = 0. The criterion for inner accuracy after m iterations has the following 
form i 

< a 



outer 

V ( t,\ II ~ " II n II ' 

II ^'outerK.'^j II L II ^outery^j II. 

where a and (3 are chosen to be 0.1 and 0.5 in our numerical simulations, since the 
rate of convergence is not sensitive to a and [3 [|, |17 



We describe next the discretization of the spatial derivatives in Eqs. ®-(|T0[) for 
the geometry of interest. Since both the Laplacian and biharmonic operators have a 
singularity at the origin of a polar coordinate system, we have found it convenient 
to use a Cartesian coordinate system and approximate the boundary conditions. We 
have used the usual 5- and 13- point discretizations of the Laplacian and biharmonic 
operators on a square grid of N x N points, which is second-order accurate in the 
mesh spacing. We approximate the circular boundary conditions on ip and C, by taking 
ip{'r,t) = ^{r,t) = for ||r|| > D/2, where is the location of a node with respect 
to the center of the domain of integration, and D is the diameter of the circular 
domain. This approximation is presumably not dangerous since boundary effects are 
known to affect the bulk behavior over length scales of the order of e~i |]TE[, which is 
large compared to the spatial discretization. Since the convective pattern has locally 
periodic solutions with wavelength 27r, except near the boundary, 8 or 12 grid points 
per wavelength are used in our numerical calculations. 

Two different kinds of initial conditions are used. The first one is a gaussian 
distributed random initial condition with zero mean value and small variance (see 
below) for ip, and ^=0. This initial condition models the conduction state. The 
second type of initial conditions used is a solution from a previous run (ipito), 'C(^o))- 
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This, for example, would correspond to studying the transition from a hexagonal state 
to a roll state in an experiment that suddenly increased the Rayleigh number. 

4 Numerical results 

In this section, we report the results of our numerical calculations based on Eqs. 

^ large aspect ratio cell, and compare the pattern evolution obtained with 
the experiments of Bodenschatz et al. We begin with the formation of a convecting 
state of hexagonal symmetry from the conducting state. Next, we present numerical 
evidence for the spontaneous formation of a rotating spiral pattern during the hexagon 
to roll transition. We also report numerical results on the formation of rotating spiral 
patterns during the transition from conduction to rolls. 

4.1 Nucleation of a pattern with hexagonal symmetry 

We consider as initial condition ■ip{r,t = 0) a Gaussian random variable with zero 
mean and variance 10^^. The forcing field chosen is /(r) = 0, simply because there 
is no influence from the lateral boundaries before the nucleated pattern reaches the 
boundary. We also neglect in this case the mean flow field. We numerically solve Eq. 
(0) in a square domain of side L = 1287r (in our dimensionless variables, this corre- 
sponds to an aspect ratio T = L/n = 128). The differential equation is discretized on 
a square grid of 512 x 512 nodes. We use g2 = 0.35 and Ax= 7r/4.25. We take e = 0.01 
except in a small square region near the center of the cell (of size 16 x 16 nodes) where 
e = 0.055. This space dependent e models a small localized inhomogeneity in one of 
the cell plates. According to the calculation by Busse and our estimate of the 
values of the parameters in the GSH equation (see the Appendix), a hexagonal pat- 
tern should be stable for both e = 0.01 and e = 0.05. The temporal evolution of the 
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pattern is shown in Fig. 1. It presents an early transient behavior during which a 
local convective region with hexagonal symmetry has just nucleated. Six fronts of 
rolls, traveling away from the hexagonal patch located at the center, propagate into 
the conduction region. This is qualitatively similar to the experimental observations 
of Bodenschatz et al. H, [l^. The observation of nucleation and growth is especially 
interesting since it provides an example of competition among different symmetries, 
i.e., a uniform conduction state as the background state, a region of hexagonal sym- 
metry being nucleated and rolls in the front region separating the two. This situation 
is also interesting from the point of view of pattern selection during front propagation 
in dimensions higher than one. It is worth pointing out, however, that the shape of 
the front region obtained numerically is somewhat different than the experimental 
one. This difference may be attributable to the fact that the numerical value e(r) 
used here is much larger than the experimental value. Unfortunately, it is very time 
consuming for us to solve the equation for the experimental value eexp ~ 10~^. 

We have also estimated the speed of propagation of the front that separates the 
hexagonal pattern and the uniform state. This speed, at the center of the planar 
sides, and along their normal direction, is approximately constant in time and equals 
v± = 0.37. The value given by marginal stability theory for the one dimensional 
Swift-Hohenberg equation ((72 = 0) is vms = 0.397 |]T9|, |2^. It will be interesting to 
measure the front velocity in the experiment. 

4.2 Conduction to hexagon transition 

In what follows, we consider a circular cell of radius r = 327r, which corresponds to 
an aspect ratio of F = r/vr = 32. A grid with A^^ nodes has been used with spacing 
Ax = Ay = 647r/N, and = 256, with an approximate boundary conditions to mimic 
a circular cell as we discussed in section 3. The initial condition ilj{r,t = 0) is again 
a random variable, gaussianly distributed with zero mean and variance 10"^ In this 
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case e = 0.1, which is also within the region in which a hexagonal pattern is stable. 
The forcing field / = everywhere except at the nodes adjacent to the boundary 
where / = 0.1. This value is of the same order as a previous estimate obtained for a 
similar experimental setup, such that the convective heat current measured during a 



ramping experiment agreed with the numerical solution of the SH equation [|l^, [16 



Figure 2 presents the evolution of the hexagonal pattern evolution from random initial 
condition. Mean fiow field effects have also been included in the calculation. 



4.3 Hexagon to roll transition, and formation of rotating spi- 
rals 

We have used the configuration shown in Fig. 2(d) as the initial condition, with 
exactly the same parameters and forcing field / as before (/ = 0.1), but have increased 
e very slowly up to e = 0.3: e = 0.1 + 1.67 x IQ-H for < t < 1200, and e = 0.3 
for t > 1200. Figure 3 shows a sequence of configurations during the early transient 
regime of the hexagon to roll transition. Rolls appear in the vicinity of the sidewall 
and tend to orient themselves parallel to it. Defects glide toward each other and invade 
nearby regions of hexagonal symmetry to create a region of rolls that spreads across 
the cell as the transition proceeds. The formation of spirals is already noticeable in 
Fig. 3(d). 

Fig. 4 shows a continuation of the sequence of configurations shown in Fig. 3. 
Fig. 4(a) shows how rolls bend rapidly to form a locally disordered texture near the 
upper left portion of the system. Defects in the disordered texture migrate away 
or annihilate each other, leaving a roughly uniform patch of rolls, and eventually 
annihilating themselves, ending in a three- armed spiral (Figs. 4(c)-(e)). The final 
state of a rotating spiral (Figs. 4(e) and (f)) is remarkably similar to the one observed 



in the experiments, and occurs at t ~ 49000 ~ 12 horizontal diffusion times [21]. The 



corresponding experimental times are in the range of 10 to 20 horizontal diffusion 
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times. 



4.4 Roll to hexagon transition 

Figures 5(a)-(d) show the evolution from a initial rotating spiral to a hexagonal 
pattern. We have used the configuration shown in Fig. 4(f) as the initial condition, 
have kept the same parameters and forcing field as before {g2 = 0.35, Qm = 50.0, 
= 10.0 and Pr=1.0), but have set e = 0.1. Initially, regions of hexagonal symmetry 
emerge near the tip and the end of the spiral, forming local domains of hexagons 
and then spreading across the cell. An interesting feature worth noticing during this 
transient period is the relative orientation between domains with different symmetries. 
If one uses the outer layer of hexagons in the domain to define a boundary, one can 
see from Fig. 5 that there is a tendency for these boundaries to be perpendicular 
to the direction defined by the rolls. In systems in which the evolution is governed 
solely by a Lyapunov functional, theoretical arguments have been given to explain 
analogous orientation phenomena Although these arguments do not apply in the 
present case, the patterns observed both in the experiments and in our numerical 
calculations still show that, locally, the boundary separating regions of hexagons and 
rolls tends to be perpendicular to the rolls. Perhaps a coupled Newell- Whitehead- 
Segel model (which has an associated Lyapunov functional) could be used to describe 



the formation of domain boundaries between hexagons and rolls |22]. 



4.5 Conduction to roll transition 



We study in this section the formation of a set of convective roll directly from the 
conducting state. We use a random initial condition (gaussianly distributed with zero 
mean and a variance of 10"^), and set e=0.3. Fi gure 6 shows that concentric rolls 
are created near the sidewall and propagate inward. Small and large length scale 
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defects anneal out rapidly leaving a disordered structure at the center of the cell (Fig. 
6(a) and (b)). Further evolution involves the annihilation of defects at the center of 
cell. Figures 6(c) and (d) show a two- armed rotating spiral in which all defects have 
been eliminated from the center of the cell. This calculation shows the importance 
of sidewall forcing and the geometric shape of the container for the formation of a 
rotating spiral. 

Another interesting feature observed in the experiments that our model can also 
reproduce is that a stable, ra-armed spiral tends toward one with fewer arms when e 
is decreased. With the same random initial condition and parameters as in Fig. 6, we 
have obtained a two-armed spiral for e = 0.3 (Fig. 6), a one-armed spiral for e = 0.26 
(Fig. 7), and a zero-armed spiral (concentric rolls) for e = 0.22. 



4.6 Stability of the rotating spiral pattern 



We discuss in this section the condition under which the rotating spiral is stable. 



Earlier work |23] established that a spiral pattern can be spontaneously formed in 
the absence of mean flow {gm = 0), but retaining the non-Boussinesq contribution 
{92 > 0). The spiral pattern obtained, however, is stationary. The addition of mean 
flow effects is sufficient to spontaneously produce a uniformly rotating spiral. We 
address here the stability of an already formed rotating spiral with respect to changes 
in the various terms of the GSH equation. 



Once the rotating spiral has been formed, we set (72 = 0.0. We use the configu- 
ration shown in Fig. 7(b) as the initial condition, with exactly the same parameters 
(e = 0.26, = 10.0, gm = 50.0, Pr = 1.0) and forcing field / as in Fig. 7(b). Figure 
8 shows the resulting evolution of the pattern. The defect in the arms of the spiral 
propagates inwards as the one-armed spiral evolves to a set of concentric rolls. This 
result suggests that non-Boussinesq effects are needed to stabilize the spiral struc- 
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ture, even under the presence of mean flow. It would be very useful to replicate this 
observation experimentally. Further, if / = rolls tend to align themselves normal 
to the sidewall, and a pattern of almost straight and parallel rolls obtains. 

In summary, sidewall forcing and non-Boussinesq effects {g2 > 0) are essential 
to produce a spiral pattern. A rotating spiral is obtained when ^ 0. Once the 
rotating spiral is formed, if g2 — 0,gm 0, the spiral decays to concentric rolls. 

4.7 Convective current versus the number of arms in a spiral 
pattern 

We have observed that for a given value of e, the number of arms of the rotating spiral 
depends on the initial condition. In order to compare the convective heat current for 
the same value of e, but for spiral patterns containing different numbers of arms, 
we proceed in the following way: the configurations shown in Fig. 7a (zero armed 
spiral for e = 0.22), Fig. 7b (one-armed spiral for e = 0.26), and Fig. 6d (two-armed 
spiral for e = 0.30) are taken as initial conditions and e is set equal to e = 0.26. 
The three spiral pattern remain stable in all cases. We then calculate the convective 
heat current as a function of the number of arms in a spiral pattern. In Fig. 9, we 
compare the spatial and temporal average of the convective current < Jconv > for 
the three cases discussed above. It is apparent from Fig. 9 that < Jconv > decreases 
with increasing the number of arms in the spiral. Since the Nusselt number is related 
to the convective current by Nu— Jconv + 1, it would also be interesting to check this 
observation experimentally. 
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5 Conclusions 



We have investigated a model of convection for non-Boussinesq fluids that allow 
patterns of various symmetries. The model used is a generalization of the Swift- 
Hohenberg equation that includes a quadratic term and coupling to large scale mean 
flows. The parameters in the equation have been chosen to match the experiments of 
Bodenschatz et al. on CO2 gas 0. An appropriate value for the control parameter 
e takes the conduction state to an ordered hexagonal state analogous to the ones 
observed in experiments. We then show that upon increasing e, the hexagonal state 
evolves into a new roll state that contains a rotating spiral pattern. The time scale 
to form a rotating spiral is on the time scale of 10 horizontal diffusion time. These 
results are also in good agreement with the experimental studies on CO2 gas. The 
observation of the stationary rotating spiral pattern is not in predicated by Busse 
1^. According to Busse when the control parameter exceeds the bifurcation point 
of the hexagon-roll state, the system is expected to evolve to a stationary parallel 
roll state. This is because they studied the convective fluid in an infinite cell. We 
have given a preliminary study of the mechanism for spiral pattern formation. Our 
calculations illustrate the strong influence of non-Boussinesq effects, sidewall forcing, 
and mean flow on the appearance and stability of the rotating spiral patterns. Our 
results for convective current of different armed-spirals show that the final state of 
the rotating spiral dependents on the initial configuration. We have seen that a zero- 
armed, a one-armed and a two-armed rotating spiral pattern could exist for the same 
control parameter. This suggests the existence of multi-states for a given control 
parameter. Our numerical results show that the two dimensional generalized Swift- 
Hohenberg quation can describe quantitatively detail the three dimensional convective 
dynamics of a fluid beyond the Boussinesq approximation. 

We conclude by suggesting some further analytical and numerical studies. It 
would be very useful to study how the instabilities, such as the zigzag, cross-roll, and 
Eckhaus instabilities are affected by the influence of the non-Boussinesq effect and the 
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finite size of tlie boundaries. It would be also very interesting to study the dynamics 
of the core and the tip (dislocation) of a rotating spiral, and to determine the speed 
of climbing motion of the tip. 
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Appendix 



In this appendix a detailed analysis of the various stationary solutions of the gener- 



alized Swift-Hohenberg equation is presented. In dimensionless units [^, the gener- 
alized SH model can be written as. 



To- 



dijj{-r, t) 
dt 



4g; 



We set 0, 



where 6j = qj ■ r, Oj = and Aj is a complex amplitude. 



(20) 



(21) 



If we substitute ipy "^^ and ip^ into Eq.(^) and keep the lowest order in the 
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amplitudes, we obtain 



To 



To 



To- 



dt 

dt 

dA3_ 
dt 






I 


OXi 


2gc 8yi 


d 


i 8^ 


8X2 


'^qc 8yl 


d 




8x3 


2gc 8yl 



Al-aA;A*-6^l(|A2^+|^3^)-c>ll|Al|^ (22) 



A2-aA\Al-hA2{\A^\''+\A3\^)-cA2\A2\\ (23) 



As-aAtA^-fc^ad^r+l^sH-c^l^r, (24) 



where a — -\/25'2, ^ — 3g'3 and c—^g^i and, 

8 ^8 . n 8 

= cosOj— + sm9j — , 
8xj 8x 8y 

8 ■ n 9 „ 8 

— = -sinOj— + cosOj—, 
8yj 8x 8y 

^1 = 0, 02 "''^'^ ^^'^y- 

If we assume a uniform solution, we simply have, 
Ml 



To 

To 

To 



8t 

8A2 
8t 

dA, 
8t 



eAi 


— aA'^A'^ 


-bA^{\A2 




-cA^\A^\^ 


eA2 


- aAlAl 


-bA2{\A, 




-cA2\A2\^ 


eAs 




-bAs{\A, 


'+\A3\') 


-cAs\As\^ 



This set of equations can be written in variational form as, 

8Ai 5C 

To- 



rn SA* 

where the Lyapunov functional C is given by. 



-e(| + + IA3P)) + aiAlAlAl + AA2A3) 



(25) 

(26) 
(27) 

(28) 

(29) 
(30) 

(31) 



(32) 
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The dynamical system Eqs. (]28|)-(pOD has three stationary and homogeneous so- 
lutions: 



Conduction state 



Ai = A2 = A3 = 0. 



Hexagonal state : 



Ai = \A\e''^\A2 = \A\e'^^,As = |A|e^^3^ 
61 + 62 + 6; = 0, 



(33) 



A\ = {-a - Ja^ + 4(26 + c)e)/(2(26 + c)). 



Roll state 



Ai = \A\e''\ A2 = /l3 = 0, 



\A\ = Je/c. 



The linear stability of these solutions is determined by the eigenvalues of the 
matrix 6'^C/6Ai6Aj, linearized around the stationary solutions. For < e < 0, both 
conduction and hexagons are stable; for < e < e^, only hexagons are stable; for 
Cr ^ e ^ both hexagon and roll are stable; and for e > e?,, only rolls are stable. 
From the corresponding values of the stability boundaries obtained for the amplitude 
equation p^, we find, 

e. = -— — = — (34) 



(86 + 4c) 15(73' 

{h-cY 3(73' 
a2(6 + 2c) 16c/| 



(6 - cf 



^93 



(35) 
(36) 



Equations 

by the Busse 
equation as 



5S| ) can be compared with the corresponding equation obtained 
This allows to determine the coefficients a, b, c in the amplitude 



3P^ 

Rr 



(37) 
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\ \ (38) 
C = ^r (39) 



Here 3?/i, 3?^ and P are the constants. In the case a rigid-rigid boundary layer, we 
have 

Rc = 1707, (40) 
= 0.8936 + 0.04959Pr"^ + 0.06787Pr"^ (41) 
= 0.69942 - 0.00472Pr-^ + 0.00832Pr-2. (42) 
where Pr is the Prandtl number, P is the non-Boussinesq parameter |Q. 

Since the system of equations (p8| - ^ ) has an associated potential £, the absolute 
stable state corresponds to the global minimum of £, while metastable states corre- 
spond to local minima. The existence of the Lyapunov functional ensures that two 
stable phases can coexist only when they have the same value of C. We define to 
be the value for which hexagons and the conduction state coexist, and ey' the value 
fot the coexistence of hexagons and rolls. We obtain, 

= -^6., (43) 

= {2b + c)\A\^ + a\A\, (44) 

where |^| is the solution of 

3c(26 + c)\A\^ + 2ac|Ap - (26 + c)\A\ - a = 0. (45) 



In order to obtain the values of the coupling coefficients g2 and g^, we calculate 
the convective current Jconv for the various patterns. Jconv is given by O], 



Jr. 



(46) 
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For a hexagonal pattern Jconv is, 



J conv 3 Ajigx 



a + Ja'^ + 4(26 + c)e /{4b + 2c) 



(x/2(^|/y^)V2 ^ /2^2/^3 + 30e)/(15V^) 



For a roll pattern Jconv is 



T _ I 4 |2 

conv \-^roU\ 



e _ 2e 
c 35(3' 



(47) 



(4^ 



By substituting the experimental values § of the threshold for ^ —2.3 x 10^^ 
in the previous expressions for the stability boundaries of the various patterns, we 
obtain, gH — 0.0345. Furthermore, at e = 0.02, we have (from Fig. 2 in [^), that 
er ^ 0.06; hence 5'|/5'3 ~ 0.045. Also, ~ 0.22, therefore, 5'|/5'3 ~ 0.0413. We have 
used g2/g3 = 0.04 in our calculations. 



By fitting the experimental convective current at e = 0.11 (from Fig. 1 in 0), we 
obtain, {Jconv)roii — 0.16 = or g^ 0.458. By fitting the experimental convective 
current at e = 0.02 (from Fig. 1 in 0) and by using 5'|/5'3 = 0.045, we obtain, 

1 2 



conv I hex 



0.04 = 3 



iV2igl/^Y/' + J2gl/g; + 30e)/(15V^) 



(49) 



or, gs ~ 0.426. 



We next discuss the calculation of the numerical parameters in the GSH equation 
related to large scale mean flow. The dynamical equations are. 



and, 



d 



ip - g2i)^ - ^3^^ 



^ - Pr(V2 - 6^) 

dt ^ ' 



where mean flow velocity f/. 



(50) 



(51) 



(52) 
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Here e = — 1 is the reduced Rayleigh number. Rc is the critical Rayleigh number 
for an infinite system, and Pr is the Prandh number. The constants tq and are 
the characteristic time and length scales, qc is the critical wave number and (72, S'a 
are the nonlinear coupling constants. 6^ is an unknown constant. Here Qm = Rc < 
Uo±{zY > /{q^ < Uoz{z)(^o{z) >). We have used the fact that near onset [ITTI , 

[V^{f,z,t),V,{f,z,t),9{f,z,t)] ^ ^[uo^{z)d^,uo,{z),9o{z)]ij{f,t) (53) 

Here V = {V±,Vz) are the velocity field and 6 is the deviation of the temperature 
from the linear conduction profile. The functions Uq±{z), Uqz{z) and 6q{z) are the first 
eigenmodes in the vertical direction of the order parameter. Here r* denotes the two- 
dimensional horizontal coordinate. The constant C = 



'<Uo,{z)eo{z)> /R,. The 
symbol < > means here an average over the vertical direction, now rescale •r' = q^f, 



2^. The rescaled GSH equation with large 



scale fiow field becomes. 



dip' 



dt 



V 



1 2 



1.1 3 



d_ 
dt' 



Pr'{V 



I 2 



J 2\ 



v'Y= v'(v'V)xvy 



(54) 
(55) 



J 2 



where g'^ = (2(72) /(^ogcv^), g'm = {^ri9mql)/ig3^i), Pr' = {ATo/^i)Pr, and c' 
h^/ql. From the experiments on CO2 g], we use Rc = 1707, qc = 3.117, = 0.148, Pr = 
l,ro = 0.07693(ro = ^in^i^p^"^ for rigid-rigid boundary conditions), gm = 2.52 (for 



rigid-rigid boundary conditions), 5'|/5'3 = 0.04, and g^ = 0.458. We finally have 
g!2 = 0.335, Pr' = 2.079, g'^ = 8.548, e' = 2.7818e, and c' ^ is an unknown constant. 



To recapitulate, we use in our numerical solution g'2 = 0.35, g'^ = 50, d 



! 2 



10 
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and Pr' = 1.0. For these values the stabihty boundaries of the various patterns are, 

e'^ = -0.008167, or, e„ = -0.00293, 

e; = 0.1633, or, = 0.0587, 

e'fe = 0.6533, or, Cf, = 0.2348, (56) 

e'y = -0.00726, or, ex = -0.00261, 

e'y, = 0.24, or, er' = 0.0863. 
Note that we have omitted the primes in the Equations (^^. 
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Figure captions 

Figure 1. Nucleation of a pattern of hexagonal symmetry in a square cell of aspect 
ratio r = 128. The values of the parameters used are g2 = 0.35, e = 0.01 and / = 0. 
In a small square region at the center of the cell, e = 0.055. The time shown ist = 611 
and dark (white) areas represent regions in which is positive (negative). 

Figure 2. Hexagonal pattern obtained from a random initial condition in a cylin- 
drical cell of aspect ratio F = 64. The values of the parameters used are g2 = 0.35, 
Qm = 50 and e = 0.1. A forcing field localized at the boundary / = 0.1 has been used. 
Four different times, (a), t = 1.3; (b), t = 51.2; (c), t = 1049.7; and, (d), t = 4649.7 
are shown. 

Figure 3. This figure shows the early stages of hexagon to roll transition produced 
by suddenly changing e from e = 0.1toe = 0.3, ina cylindrical cell of aspect ratio 
r = 64. The initial condition is the uniform hexagonal pattern shown in Fig. 2a. 
Four different times, (a), t = 480; (b), t = 720; (c), ^ = 840; and, (d), t = 960 are 
shown. Rolls appear near defects and sidewall boundaries and spread through the 
cell as the transition proceeds. 

Figure 4. Formation of a rotating three-armed spiral in a cylindrical cell of aspect 
ratio r = 64, with g2 = 0.35, gm = 50, e = 0.3 and / = 0.1. The configurations 
shown are continuation of those in Fig. 3. The times shown are (a), t=1400; (b), t= 
7440; (c), t=18000; (d), t=41040; (e), t=48240; and, (f), t=64080. The final rotating 
spiral pattern is shown in (e) and (f) is stable. 

Figure 5. Formation of a hexagonal pattern from a three-armed rotating spiral 
(Fig. 4(f)), with g2 = 0.35, g^ = 50, e = 0.1 and f=0.1. An interesting feature to 
note is that the orientation of the hexagonal domains tends to be normal to the rolls. 
Four configurations during th early transition period are shown here: (a), t=649; (b), 
t=1299; (c), t=2130; and, (d), t=3090. 
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Figure 6. Formation of a two-armed rotating spiral from random initial conditions, 
with §2 = 0.35, Qm = 50, e = 0.3 and / = 0.1. The configurations are shown at (a), 
t=70.1; (b), t=190.7; (c), t=310.7; and, (d), t=1990.7. 

Figure 7. Formation of a zero-armed, and a one-armed rotating spiral pattern 
obtained from a random initial condition with different values of e from Fig. 6. All 
the other parameters are the same as in Fig. 6. (a), a zero-armed spiral with e = 0.22, 
and (b), a one-armed spiral with e = 0.26. 

Figure 8. Study of the stability of the spiral pattern. Wc use the one-armed spiral 
shown in Fig. (7b) as initial condition, switch g2 from gf2=0.35 to g2=0.0 and keep 
all the other parameters unchanged. Concentric rolls propagate inwards as the core 
of the one-armed spiral starts shrinking. Finally the one-armed spiral decays to a set 
of concentric rolls. The times shown are (a), t=0; (b), t= 180; (c), t= 252; and, (d), 
t= 288. 

Figure 9. Convective heat current Jconv versus the number of arms for the same 
value of e = 0.26 The average convective current decreases with increasing the number 
of arms in the spiral. 
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